ppseq: An R Package for Sequential Predictive Probability Monitoring

Introduction

In the context of cytotoxic treatments, phase 1 trials in oncology traditionally have the primary aim of identifying the maximum tolerated dose (MTD), defined as the highest dose that still maintains a certain pre-specified rate of toxicity, most commonly 30%, in a dose-escalation phase. Designs for dose-escalation trials include the rule-based 3+3 design and the model-based continual reassessment method, among others. But with increasing study focused on non-cytotoxic treatments, such as immunotherapies, the MTD either may not exist or may not be relevant, and toxicities may either develop much later or even be chronic, making these treatments difficult to study with traditional dose-escalation designs (Pestana et al. 2020). As a result, it is becoming increasingly common to include dose-expansion cohorts, in which additional patients are enrolled in phase 1 after the dose-escalation phase is complete. In this setup, the dose-escalation phase is considered phase 1a and used to assess the initial safety of multiple doses, then the dose-expansion phase is considered phase 1b and can have a variety of aims including to further refine the safety of one or more doses, to assess preliminary efficacy, to explore the treatment in various disease-specific subtypes that all share a common biomarker that the treatment is targeting, or to further characterize the pharmacokinetics and/or pharmacodynamics. The use of dose-expansion cohorts increased from 12% in 2006 to 38% in 2011 (Manji et al. 2013) and trials with dose-expansion cohorts led to higher response rates and more frequent success in phase 2 trials (Bugano et al. 2017).

But dose-expansion cohorts are not always planned in advance, and therefore may be subject to on-the-fly decision making that can lead to large sample sizes and poor statistical properties. For example, the KEYNOTE-001 trial of pembrolizumab, initially designed as a 3+3 dose-escalation trial, included multiple protocol amendments and ultimately enrolled a total of 655 patients across five melanoma expansion cohorts and 550 patients across four non-small-cell lung cancer expansion cohorts (Khoja et al. 2015). In a basket trial of atezolizumab, an anti-PD-L1 treatment in patients with a variety of cancers and both with and without PD-L1 expression, an expansion cohort in metastatic urothelial bladder cancer ultimately enrolled 97 patients and evaluated 95, despite the fact that no expansion cohort in this disease subtype was originally planned in the trial protocol. The expansion cohort in metastatic urothelial carcinoma was rather added later in a protocol amendment in which the sample size was increased from what was initially planned (Petrylak et al. 2018; Powles et al. 2014).

Bayesian sequential predictive probability monitoring provides a natural framework for early phase oncology trial design, allowing flexibility to stop early for futility or safety concerns while also maintaining traditional levels of power and control of type I error. As a result, predictive probability has been proposed as an approach to interim monitoring in clinical trials previously (Dmitrienko and Wang 2006; Lee and Liu 2008; Saville et al. 2014; Hobbs, Chen, and Lee 2018). However, in order to be useful to investigators designing such trials, software must be made available to select an optimal design. This paper introduces the ppseq package for the R software language (R Core Team 2020), which provides functions to design one-sample and two-sample early phase clinical trials using sequential predictive probability monitoring for futility. Interactive plots produced using the ggplot2 package (Wickham 2016) and the plotly package (Sievert 2020) compare designs based on different thresholds for decision making. Moreover, we demonstrate optimization criteria for selecting an ideal predictive probability monitoring design using the ppseq package. While the ppseq package was developed with early phase oncology clinical trials in mind, the methodology is general and can be applied to any application of sequential futility monitoring in clinical trial design.

Predictive probability monitoring

Consider the setting of a binary outcome, such as tumor response as measured by the RECIST criteria in the setting of a study in patients with solid tumors, where each patient, denoted by \(i\), enrolled in the trial either has a response such that \(x_i = 1\) or does not have a response such that \(x_i = 0\). Then \(X = \sum_{i=1}^n x_i\) represents the number of responses out of \(n\) currently observed patients up to a maximum of \(N\) total patients. Let \(p\) represent the probability of response, where \(p_0\) represents the null response rate under no treatment or the standard of care treatment and \(p_1\) represents the alternative response rate under the experimental treatment. Most dose-expansion studies with an efficacy aim will wish to test the null hypothesis \(H_0: p \leq p_0\) versus the alternative hypothesis \(H_1: p \geq p_1\).

The Bayesian paradigm of statistics is founded on Bayes rule, which is a mathematical theory that specifies how to combine the prior distributions that define prior beliefs about parameters, such as the response rate \(p\), with the observed data, such as the total number of responses \(X\), yielding a posterior distribution. Here the prior distribution of the response rate \(\pi(p)\) has a beta distribution \(Beta(a_0, b_0)\) and our data \(X\) have a binomial distribution \(bin(n, p)\). Combining the likelihood function for the observed data \(L_x(p) \propto p^x (1-p)^{n-x}\) with the prior, we obtain the posterior distribution of the response rate, which follows the beta distribution \(p|x \sim Beta(a_0 + x, b_0 + n - x)\). A posterior probability threshold \(\theta\) would be pre-specified during the trial design stage. At the end of the trial if the posterior probability exceeded the pre-specified threshold, i.e. if \(\Pr(p>p_0 | X) > \theta\), the trial would be declared a success.

The posterior predictive distribution of the number of future responses \(X^*\) in the remaining \(n^*=N-n\) future patients follows a beta-binomial distribution \(Beta-binomial(n^*, a_0 + x, b_0 + n - x)\). Then the posterior predictive probability (PPP), is calculated as \(PPP = \sum_{{x^*}=0}^{n^*} \Pr(X^*=x^*|x) \times I(\Pr(p>p_0 | X, X^*=x^*) > \theta)\). The posterior predictive probability represents the probability that, at any given interim monitoring point, the treatment will be declared efficacious at the end of the trial when full enrollment is reached, conditional on the currently observed data and the specified priors. We would stop the trial early for futility if the posterior predictive probability dropped below a pre-specified threshold \(\theta^*\), i.e. \(PPP<\theta^*\). Predictive probability thresholds closer to 0 lead to less frequent stopping for futility whereas thresholds near 1 lead to frequent stopping unless there is almost certain probability of success. Predictive probability provides an intuitive interim monitoring strategy for clinical trials that tells the investigator what the chances are of declaring the treatment efficacious at the end of the trial if we were to continue enrolling to the maximum planned sample size, based on the data observed in the trial to date.

Package overview

The ppseq package facilitates the design of clinical trials utilizing sequential predictive probability monitoring for futility. The goal is to establish a set of decision rules at the trial planning phase that would be used for interim monitoring during the course of the trial. The main challenge in designing such a trial is joint calibration of the posterior probability and posterior predictive probability thresholds to be used in the trial in order to maintain the desired level of power and type I error. The function will evaluate a grid of posterior and predictive thresholds provided by the user as vector inputs through the argument for posterior thresholds and the argument for predictive thresholds. Other required arguments include the null response rate , the alternative response rate , a vector of sample sizes at which interim analyses are to be performed as well as the maximum total sample size . Default arguments are available for the direction of the alternative hypothesis , a vector of the two hyperparameters of the prior beta distribution , and the number of posterior samples and the number of simulated trial datasets , but customization is also possible. The additional argument can specify the clinically meaningful difference between groups in the case of a two-sample trial design. The function returns a list, the first element of which is a containing the posterior threshold, predictive threshold, the mean sample size under the null and the alternative, the proportion of positive trials under the null and alternative, and the proportion of trials stopped early under the null and alternative. The proportion of positive trials under the null represents the type I error and the proportion of positive trials under the alternative represents the power. The option will print the results summary for each combination of thresholds, filtered by an acceptable range of type I error and minimum power, if desired.

Optimization

After obtaining results for all combinations of evaluated posterior and predictive thresholds, the next step is to select the ideal design from among the various options. The ppseq package introduces two optimization criteria to assist users in making a selection. The first, called the “optimal accuracy” design, identifies the design that minimizes the Euclidean distance to the top left point on a plot of the type I error by the power. The second, called the “optimal efficiency” design, identifies the design that minimizes the Euclidean distance to the top left point on a plot of the average sample size under the null by the average sample size under the alternative, subject to constraints on the type I error and power. The function will return a list that contains the details of each of the two optimal designs.

Decision rules

To ease the implementation of clinical trials designed with sequential predictive probability monitoring, once a design has been selected, a table of decision rules can be produced using the function. The function takes the sample sizes at which interim analyses are to be performed as well as the maximum total sample size , the null value to compare to in the one-sample case (set to in the two-sample case), the posterior threshold of the selected design , and the predictive threshold of the selected design . Default arguments are provided for the direction of the alternative hypothesis , a vector of the two hyperparameters of the prior beta distribution , and the number of posterior samples , though customization is also possible. The additional argument can specify the clinically meaningful difference between groups in the case of a two-sample trial design (set to in the one-sample case). The function results in a . In the one-sample case, the includes a column for the sample size at each interim analysis, a column for the decision point \(r\), and a column for the posterior predictive probability associated with that decision . The trial would stop at a given look if the number of observed responses is less than or equal to \(r\) otherwise the trial would continue enrolling if the number of observed responses is greater than \(r\). At the end of the trial when the maximum planned sample size is reached, the treatment would be considered promising if the number of observed responses is greater than \(r\). In the two-sample case, the includes columns for and , the sample size at each interim analysis in the control and experimental arms, respectively. There are also columns for and , the number of responses in the control arm and experimental arm leading to the decision, respectively. Finally, there is a column for the posterior predictive probability associated with that decision .

Visualizations

Finally, to assist users in comparing the results of the various design options, a option is available for the results of that allows creation of static plots using the ggplot2 package (Wickham 2016) or interactive plots using the plotly package (Sievert 2020). Two plots are produced, one plotting type I error by power and indicating the optimal accuracy design, and one plotting the average sample size under the null by the average sample size under the alternative and indicating the optimal efficiency design. The motivation for including an interactive graphics option was the utility of the additional information available when hovering over each point. Instead of simply eyeballing where points fall along the axes, users can see the specific type I error, power, average sample size under the null, average sample size under the alternative, the posterior and predictive thresholds associated with the design, as well as the distance to the upper left point on the plot. A option is also available for the results of . In the one-sample case it produces a single plot showing the sample size at each interim analysis on the x-axis and the possible number of responses at each interim analysis on the y-axis. In the two-sample case a grid of plots is produced, with one plot for each interim analysis. The x-axis shows the number of possible responses in the control group and the y-axis shows the number of possible responses in the experimental group. In both cases, the boxes are colored green for a “proceed” decision and red for a “stop” decision for each combination and the hover box produced by plotly provides the details.

Case study

To demonstrate the functionality of the ppseq package, we focus on a re-design of a dose-expansion cohort for the study of atezolizumab in metastatic urothelial carcinoma patients (mUC) using sequential predictive probability monitoring. Atezolizumab is an anti-PD-L1 treatment that was originally tested in the phase 1 setting in a basket trial across a variety of cancer sites harboring PD-L1 mutations. The atezolizumab expansion study in mUC had the primary aim of further evaluating safety, pharmacodynamics and pharmacokinetics and therefore was not designed to meet any specific criteria for type I error or power. An expansion cohort in mUC was not part of the original protocol design, but was rather added later through a protocol amendment. The expansion cohort in mUC ultimately evaluated a total of 95 participants (Powles et al. 2014). Other expansion cohorts that were included in the original protocol, including in renal-cell carcinoma, non-small-cell lung cancer, and melanoma, were planned to have a sample size of 40. These pre-planned expansion cohorts were designed with a single interim analysis for futility that would stop the trial if 0 responses were seen in the first 14 patients enrolled. According to the trial protocol, this futility rule is associated with at most a 4.4% chance of observing no responses in 14 patients if the true response rate is 20% or higher. The protocol also states the widths of the 90% confidence intervals for a sample size of 40 if the observed response rate is 30%. There was no stated decision rule for efficacy since efficacy was not an explicit aim of the expansion cohorts. In the re-design we assume a null, or unacceptable, response rate of 0.1 and an alternative, or acceptable, response rate of 0.2. We plan a study with up to a total of 95 participants. In our sequential predictive probability design we will check for futility after every 5 patients are enrolled. We consider posterior thresholds of 0, 0.7, 0.74, 0.78, 0.82, 0.86, 0.9, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99, 0.999, 0.9999, 0.99999, and 1, and predictive thresholds of 0.05, 0.1, 0.15, and 0.2.

First we install and load the ppseq package.

install.packages("ppseq")

library(ppseq)

We use the function to obtain the operating characteristics of designs based on on each combination of posterior and predictive thresholds. Because of the inherent computation intensity in these calculations, this function relies on the future (Bengtsson 2020) and furrr (Vaughan and Dancho 2021) packages to parallelize computations. The user will be responsible for setting up a call to that is appropriate to their operating environment and simulation setting. The example in this case study was run on a Unix server with 192 cores, and we wished to use 76 cores to accommodate the 76 distinct designs that result from the 19 posterior by 4 predictive threshold grid. Because the code takes some time to run, the results of the below example code are available as a dataset called included in the ppseq package.

library(future)

set.seed(123)

plan(multicore(workers = 76))

one_sample_cal_tbl <- 
  calibrate_thresholds(
    p_null = 0.1, 
    p_alt = 0.2, 
    n = seq(5, 95, 5),
    N = 95, 
    pp_threshold = c(0, 0.7, 0.74, 0.78, 0.82, 0.86, 0.9, 0.92, 0.93, 0.94, 
                     0.95, 0.96, 0.97, 0.98, 0.99, 0.999, 0.9999, 0.99999, 1),
    ppp_threshold = seq(0.05, 0.2, 0.05),
    direction = "greater", 
    delta = NULL, 
    prior = c(0.5, 0.5), 
    S = 5000, 
    nsim = 1000
  )

Next we print the results table using the option, and limited to designs with type I error between 0.05 and 0.1, and a minimum power of 0.7. We find that 35 of the 76 designs meet these criteria for type I error and power.

print(
  one_sample_cal_tbl,
  type1_range = c(0.05, 0.1),
  minimum_power = 0.7
)

# A tibble: 14 x 8
   pp_threshold ppp_threshold mean_n1_null prop_pos_null
          <dbl>         <dbl>        <dbl>         <dbl>
 1         0.82          0.1          42.9         0.096
 2         0.82          0.15         39.1         0.081
 3         0.82          0.2          35.3         0.076
 4         0.86          0.1          42.8         0.096
 5         0.86          0.15         39.1         0.081
 6         0.86          0.2          35.0         0.073
 7         0.9           0.05         51.1         0.072
 8         0.9           0.1          38.6         0.061
 9         0.9           0.15         35.0         0.057
10         0.92          0.05         50.8         0.071
11         0.92          0.1          38.6         0.061
12         0.92          0.15         35.1         0.058
13         0.93          0.05         50.1         0.05 
14         0.95          0.05         46.2         0.05 
# ... with 4 more variables: prop_stopped_null <dbl>,
#   mean_n1_alt <dbl>, prop_pos_alt <dbl>, prop_stopped_alt <dbl>

We use the function to obtain the details of the optimal accuracy and optimal efficiency designs, limited to our desired range of type I error and minimum power. We find that the optimal accuracy design is the one with posterior threshold 0.9 and predictive threshold 0.05. It has a type I error of 0.072, power of 0.883, average sample size under the null of 51, and average sample size under the alternative of 91. The optimal efficiency design is the one with posterior threshold of 0.92 and predictive threshold of 0.1. It has a type I error of 0.06, power of 0.796, average sample size under the null of 39, and average sample size under the alternative of 82. For comparison, the original design of the atezolizumab expansion cohort in mUC, with a single look for futility after the first 14 patients, has a type I error of 0.005, power of 0.528, average sample size under the null of 76, and average sample size under the alternative of 92.

optimize_design(
  one_sample_cal_tbl, 
  type1_range = c(0.05, 0.1), 
  minimum_power = 0.7
)

$`Optimal accuracy design:`
# A tibble: 1 x 6
  pp_threshold ppp_threshold `Type I error` Power
         <dbl>         <dbl>          <dbl> <dbl>
1          0.9          0.05          0.072 0.882
  `Average N under the null` `Average N under the alternative`
                       <dbl>                             <dbl>
1                       51.1                              90.2

$`Optimal efficiency design:`
# A tibble: 1 x 6
  pp_threshold ppp_threshold `Type I error` Power
         <dbl>         <dbl>          <dbl> <dbl>
1         0.92           0.1          0.061 0.796
  `Average N under the null` `Average N under the alternative`
                       <dbl>                             <dbl>
1                       38.6                              81.6

To compare these optimal designs with all other designs, we can use the function with the option to obtain interactive visualizations.

plot(
  one_sample_cal_tbl, 
  type1_range = c(0.05, 0.1), 
  minimum_power = 0.7,
  plotly = TRUE
)

Figure 1: Plot of design options made with plotly. The color represents the Euclidean distance to the top left point and the optimal design is indicated by a diamond.

Figure 1: Plot of design options made with plotly. The color represents the Euclidean distance to the top left point and the optimal design is indicated by a diamond.

In this case we may choose to use the optimal efficiency design, which has the desirable trait of a very small average sample size under the null of just 39 patients, while still maintaining reasonable type I error of 0.06 and power of 0.796. This design would allow us to stop early if the treatment were inefficacious, thus preserving valuable financial resources for use in studying more promising treatments and preventing our human subjects from continuing an ineffective treatment. Finally, we generate the decision table associated with the selected design for use in making decision at each interim analysis during the conduct of the trial. Because of the computational time involved, the results of the below example code are available as a dataset called included in the ppseq package. In the results table, we see that at the first interim futility look after just 5 patients, we would not stop the trial. After the first 10 patients we would stop the trial if there were 0 responses, and so on. At the end of the trial when all 95 patients have accrued, we would declare the treatment promising of further study if there were greater than or equal to 14 responses.

set.seed(123)

one_sample_decision_tbl <- 
  calc_decision_rules(
    n = seq(5, 95, 5),
    N = 95, 
    theta = 0.92, 
    ppp = 0.1, 
    p0 = 0.1, 
    direction = "greater",
    delta = NULL, 
    prior = c(0.5, 0.5), 
    S = 5000
  )

one_sample_decision_tbl

# A tibble: 19 x 3
       n     r     ppp
   <dbl> <int>   <dbl>
 1     5    NA NA     
 2    10     0  0.0634
 3    15     0  0.022 
 4    20     1  0.0844
 5    25     1  0.0326
 6    30     2  0.0702
 7    35     2  0.0288
 8    40     3  0.0536
 9    45     4  0.0708
10    50     4  0.0344
11    55     5  0.0478
12    60     6  0.0622
13    65     7  0.0888
14    70     8  0.095 
15    75     8  0.033 
16    80     9  0.0326
17    85    10  0.0346
18    90    11  0.0162
19    95    13  0     

plot(one_sample_decision_tbl)

Figure 2: Plot of decision rules made with plotly. The color indicates whether the trial should stop or proceed for a given number of responses at each interim analysis.

Summary

With the focus of early stage clinical trial research in oncology shifting away from the study of cytotoxic treatments and toward immunotherapies and other non-cytotoxic treatments, new approaches to clinical trial design are needed that move beyond the traditional search for the maximum tolerated dose (Hobbs et al. 2019). Bayesian sequential predictive probability monitoring provides a natural and flexible way to expand the number of patients studied in phase 1 or to design phase 2 trials that allow for efficient early stopping for futility while maintaining control of type I error and power. The ppseq package implements functionality to evaluate a range of posterior and predictive thresholds for a given study design and identify the optimal design based on accuracy (i.e. type I error and power) or efficiency (i.e. average sample sizes under the null and alternative). Interactive visualization options are provided to ease comparison of the resulting design options. Once an ideal design is selected, a table of decision rules can be obtained to make trial conduct simple and straightforward.